############################################################
## Replication code for:                                  ##
## Yuri M. Zhukov & Brandon M. Stewart                    ##
## "Choosing Your Neighbors"                              ##
## Updated Sep 2011                                       ##
## Forthcoming in International Studies Quarterly         ##
##                                                        ##
## Contains replication code for all Figures and Tables   ##
## except Figure 2 (which was generated in TeX)           ##
## ## All plots print to PDF files                        ##
############################################################

## Clear memory
rm(list=ls())

## Set working directory (where you unzipped the archive)
setwd("/MyDirectory")

## Load data and source functions
load("Data.RData")
source("Functions.R")

## Load libraries
packages(reshape)
packages(MASS)
packages(sna)
packages(foreign)
packages(maptools)
packages(Zelig)
packages(mvtnorm)
packages(mnormt)
packages(spdep)
packages(RColorBrewer) 
packages(classInt) 
packages(numDeriv)
packages(zipfR)
packages(sp)
packages(xtable)





#########################################################
#########################################################
## Finding the “True” Network: A Monte Carlo Study     ##
#########################################################
#########################################################


#######################
## Generate 1000 random graphs based on contiguity
#######################

Ws <- array(NA,dim=c(1000,159,159))
x <- seq(0,1,length.out=1000)
probs <- 1-Rbeta.inv(x, a=.25, b=.25)
for(i in 1:1000){
W <- matrix(mat_cont,159,159)
W. <- W
prob <- probs[i]
W[which(W.==1)] <- prob
W[which(W.==0)] <- 1-prob
Ws[i,,] <- W
}
RW <- rgraph(159,m=1000,tprob=Ws,mode="graph",diag=T)

## Graph Correlation matrix
CorRW <- gcor(RW)

## Save output
#save(CorRW,file="BigGCor.RData")
#load("BigGCor.RData")

#######################
## Generate fake data: SAR
#######################

## Extract W matrix
W1 <- RW[1,,]

## Row-normalize
W1 <- row.norm(W1)

## Set parameters
rho <- .1
alpha <- .01
beta <- 4
Betas <- c(alpha,beta)

## Synthetic predictors
set.seed(123)
X <- matrix(NA,nrow=159,ncol=2)
X[,1] <- 1
X[,2] <- rnorm(159)
iN <- rep(1,159)

## Error term 
set.seed(123)
u <- rmnorm(n = 1, mean = 0, varcov = diag(159)) 
u <- t(u)

## Calculate Y through autoregressive process (y = pWy + XB + u)
y <- invIrW(mat2listw(W1), rho) %*% X%*%Betas + invIrW(mat2listw(W1), rho) %*% u
rbinom(n=159,size=1,prob=pnorm(y))

## Check
ck <- cbind(y,rho*W1%*%y + X%*%Betas + u) 
if(mean(round(ck[,1],12)==round(ck[,2],12))==1){print("OK")}

## Dataset 
data. <- cbind(y,X)
data. <- as.data.frame(data.)
colnames(data.) <- c("y","X0","X1")


#######################
## MLE loop: summarize results of 1000 SAR regressions
#######################

ords <- order(CorRW[1,])
SARs <- matrix(NA,nrow=1000,ncol=6)
SARs <- as.data.frame(SARs)
names(SARs) <- c("GCor","Rho","RhoSE","AIC","BIC","SSR")
SARs$GCor <- CorRW[1,ords]
for(i in 1:1000){
	print(i)
	W <- RW[ords[i],,]
	W <- row.norm(W)
mod <- lagsarlm(y ~ X1, data = data., listw=mat2listw(W), na.action="na.omit", zero.policy=T)
SARs[i,2] <- mod$rho
SARs[i,3] <- mod$rho.se
SARs[i,4] <- -2*mod$LL + 2*mod$parameters
SARs[i,5] <- -2*mod$LL + mod$parameters*log(length(mod$residuals))
SARs[i,6] <- sum(mod$residuals^2)
}
SARs


#######################
## Cross-validation
#######################


rmse.mat <- matrix(NA,nrow=1000,ncol=100)
for(j in 1:1000){
    print(j)
W <- RW[ords[j],,]
W <- row.norm(W)
data.$Wy <- W%*%y
for(i in 1:100){
out. <- sample(1:nrow(data.),size=round(.1*nrow(data.)))
out.samp <- data.[out.,]
in.samp <- data.[-out.,]
mod <- zelig(y ~ Wy + X1, data = in.samp, model="normal",cite=F)
mod.x <- setx(mod, fn = NULL,data=out.samp)
mod.y <- sim(mod,mod.x)
mod.fit <- apply(mod.y$qi$ev,2,mean)
rmse.mat[j,i] <- sqrt(mean((out.samp$y - mod.fit)^2))
}}
rmses <- apply(rmse.mat,1,mean)
SARs$RMSE <- apply(rmse.mat,1,mean)
SARs$RMSE.l <- apply(rmse.mat,1,function(x){quantile(x,.025)})
SARs$RMSE.u <- apply(rmse.mat,1,function(x){quantile(x,.975)})


## Save output
#save(SARs,file="SARS.RData")
#load("SARS.RData")


#######################
## FIGURE 1: Monte Carlo results
#######################

## 1a: Strength of spatial autocorrelation
pdf("Fig1a.pdf",width=4,height=3)
par(mar=c(4,4,.5,.5))
plot(SARs$GCor,SARs$Rho,cex=.5,col="white",xlab="Graph Correlation",ylab=expression(hat(rho)))
segments(x0=SARs$GCor,x1=SARs$GCor,y0=SARs$Rho+1.96*SARs$RhoSE,y1=SARs$Rho-1.96*SARs$RhoSE,col="grey",lwd=2)
points(SARs$GCor,SARs$Rho,cex=.3,col="black")
abline(h=.1,col="red")	# Color
#abline(h=.1,col="black")	# B&W
abline(h=0,lty="dashed")
legend(x="bottomright",bty="n",col="red",lwd=1,legend="True value",cex=1)			# Color
#legend(x="bottomright",bty="n",col="black",lwd=1,legend="True value",cex=1)			# B&W
dev.off()

## 1b: Goodness of fit (AIC)
pdf("Fig1b.pdf",width=4,height=3)
par(mar=c(4,4,.5,.5))
plot(SARs$GCor,SARs$AIC,cex=.5,xlab="Graph Correlation",ylab="AIC",ylim=c(-1650,5))
dev.off()

## 1b (alt): Goodness of fit (BIC)
pdf("Fig1bALT.pdf",width=4,height=3)
par(mar=c(4,4,.5,.5))
plot(SARs$GCor,SARs$BIC,cex=.5,xlab="Graph Correlation",ylab="BIC",ylim=c(-1650,5))
dev.off()

## 1c: Cross-validation
pdf("Fig1c.pdf",width=4,height=3)
par(mar=c(4,4,.5,.5))
plot(SARs$GCor,SARs$RMSE,cex=.5,col="white",xlab="Graph Correlation",ylab="RMSE (out of sample)",ylim=c(0,.35))
segments(x0=SARs$GCor,x1=SARs$GCor,y0=SARs$RMSE.l,y1=SARs$RMSE.u,col="grey",lwd=1)
points(SARs$GCor,SARs$RMSE,cex=.3,col="black")
dev.off()


#########################################################
#########################################################
## Application: The Spread of Democracy                ##
## Step 1: Enumeration of Network Specifications       ##
#########################################################
#########################################################


## Create array of candidate matrices
labels <- c("Geographic (CONT)","Geographic (MDN)","Geographic (KNN4)","Geographic (SOI)","Ethnic (MDN)","Ethnic (KNN4)","Ethnic (pSOI)","Trade (MDN)","Trade (KNN4)","Trade (pSOI)","IGO (MDN)","IGO (KNN4)","IGO (pSOI)","Alliance Ties")
mat_array <- array(NA,dim=c(14,159,159))
mat_array[1,,] <- mat_cont
mat_array[2,,] <- mat_dist
mat_array[3,,] <- mat_knn4
mat_array[4,,] <- mat_soi
mat_array[5,,] <- mat_ethmd
mat_array[6,,] <- mat_ethk4
mat_array[7,,] <- mat_ethsoi
mat_array[8,,] <- mat_trmd
mat_array[9,,] <- mat_trk4
mat_array[10,,] <- mat_trsoi
mat_array[11,,] <- mat_igomd
mat_array[12,,] <- mat_igok4
mat_array[13,,] <- mat_igosoi
mat_array[14,,] <- mat_ally
mat.cor <- gcor(mat_array)
ords <- order(mat.cor[1,])



##########################
## FIGURE 3: Plot the connections
## NOTE: This code generates each connectivity plot separately.
## Collage of the 14 plots was created in TeX.
##########################

map_crd <- cbind(subdata$coor.x,subdata$coor.y)
library(maps)
data(worldMapEnv)
trs <- rep(.5,14)
trs[c(2,5,8,11)] <- .1
trs[14] <- .3

## Generate plot for each connectivity
for(j in 1:14){
pdf(paste("Fig3_",j,".pdf",sep=""),height=2,width=3)
par(mar=c(0,0,0,0))
W_mat <- mat2listw(mat_array[j,,],row.names=rownames(mat_cont))
map('world', fill = F,col="lightgray",xlim=c(-141,180),ylim=c(-55,78))
plot(W_mat,coords=map_crd,pch=19, cex=0.1, col=rgb(0,.7,.1, alpha=trs[j]),add=T,lwd=.8)
title(paste(labels[j],"\n"))
dev.off()
}


pdf("Fig3_Blank.pdf",height=2,width=3)
par(mar=c(0,0,0,0))
plot.new()
dev.off()



#######################
## FIGURE 4: Correlation Matrix (b&w)
#######################

rownames(mat.cor) <- labels
colnames(mat.cor) <- rownames(mat.cor) 

# Color palette
pal <- c("gray0","gray25","gray50","gray75","white","gray75","gray50","Gray25","gray0")
cols.mat <- matrix(NA,nrow=nrow(mat.cor),ncol=ncol(mat.cor))
for(i in 1:ncol(mat.cor)){
classes_fx <- classIntervals(var=mat.cor[,i], n=10, style="fixed",
fixedBreaks=c(-1,-.8,-.6,-.4,-.2,0,.2,.4,.6,.8,1),  rtimes = 3)
cols <- findColours(classes_fx, pal)
cols.mat[,i] <- cols[1:length(cols)]
}
cols.mat <- cols.mat[ nrow(cols.mat):1,]

# Point types
pchs <- matrix(NA,nrow=nrow(mat.cor),ncol=ncol(mat.cor))
pchs[mat.cor < 0] <- 4
pchs[mat.cor >= 0] <- NA
pchs <-pchs[ nrow(pchs):1,]
pch.leg <- c(rep(4,5),rep(NA,5))

# Plot it and save
pdf("Fig4.pdf",height=4,width=5.5)
par(mar=c(0,0,0,0))
plot(c(-5,nrow(mat.cor)+3), c(-5,nrow(mat.cor)+1), type = "n", xlab="", ylab="", asp = 1,axes=F)
for(i in 1:nrow(cols.mat)){
for(j in 1:ncol(cols.mat)){
points(x=i,y=j,cex=2,col=cols.mat[j,i],pch=15)
points(x=i,y=j,cex=2,col="black",pch=pchs[j,i])
}}
text(x=0.6,y=ncol(cols.mat):1,pos=2,labels = colnames(mat.cor),cex=.7)
text(x=2:(ncol(cols.mat)+1)-.55,y=0,pos=2,labels = colnames(mat.cor),cex=.7, srt = 90)
legend(x="right",cex=.8,col=attr(cols,"palette"),bty="n",legend=names(attr(cols, "table")),title="Correlation",ncol=1,pch=15,pt.cex=1.5)
legend(x="right",cex=.8,col="black",bty="n",legend=names(attr(cols, "table")),title="Correlation",ncol=1,pch=pch.leg,text.col=NA,title.col=NA,pt.cex=1.5)
dev.off()



#######################
## FIGURE 4 (alt): Correlation Matrix (color)
#######################

rownames(mat.cor) <- labels
colnames(mat.cor) <- rownames(mat.cor) 

# Color palette
pal <- c("darkblue","blue","lightblue3","lightblue","beige","orange","darkorange","red","darkred")
cols.mat <- matrix(NA,nrow=nrow(mat.cor),ncol=ncol(mat.cor))
for(i in 1:ncol(mat.cor)){
classes_fx <- classIntervals(var=mat.cor[,i], n=10, style="fixed",
fixedBreaks=c(-1,-.8,-.6,-.4,-.2,0,.2,.4,.6,.8,1),  rtimes = 3)
cols <- findColours(classes_fx, pal)
cols.mat[,i] <- cols[1:length(cols)]
}
cols.mat <- cols.mat[ nrow(cols.mat):1,]

# Plot it and save
pdf("Fig4_color.pdf",height=4,width=5.5)
par(mar=c(0,0,0,0))
plot(c(-5,nrow(mat.cor)+3), c(-5,nrow(mat.cor)+1), type = "n", xlab="", ylab="", asp = 1,axes=F)
for(i in 1:nrow(cols.mat)){
for(j in 1:ncol(cols.mat)){
points(x=i,y=j,cex=2,col=cols.mat[j,i],pch=15)
}}
text(x=0.6,y=ncol(cols.mat):1,pos=2,labels = colnames(mat.cor),cex=.7)
text(x=2:(ncol(cols.mat)+1)-.55,y=0,pos=2,labels = colnames(mat.cor),cex=.7, srt = 90)
legend(x="right",cex=.8,fill=attr(cols,"palette"),bty="n",legend=names(attr(cols, "table")),title="Correlation",ncol=1)
dev.off()








##########################################################
##########################################################
## Application: The Spread of Democracy                 ##
## Step 2: Structural Comparisons and Model Diagnostics ##
##########################################################
##########################################################


#############################
## Data management
#############################

## Reduce size of dataset to bare minimum needed for models
data. <- data[,c("baut","bautlag","laggdppc","energy2","pnbdem","cwar","ipyears","egr","propdem","nbtd","demtr","Wpropdem.dist.2","Wpropdem.k4.2","Wpropdem.soi.2","Wdemtr.dist.2","Wdemtr.k4.2","Wdemtr.soi.2","W.pdem.mdn.eth","W.demtr.mdn.eth","W.pdem.k4.eth","W.demtr.k4.eth","W.pdem.soi.eth","W.demtr.soi.eth","W.pdem.mdn.tr","W.demtr.mdn.tr","W.pdem.k4.tr","W.demtr.k4.tr","W.pdem.soi.tr","W.demtr.soi.tr","W.pdem.al","W.demtr.al","W.pdem.mdn.igo","W.demtr.mdn.igo","W.pdem.k4.igo","W.demtr.k4.igo","W.pdem.soi.igo","W.demtr.soi.igo")]

## Drop missing values
data. <- na.omit(data.)


###############################
## FIGURE 5a: SD shifts
###############################

## Run GW2006 model on 14 neighbor matrices
modz.sd <- models(data.)

## Calculate predicted shifts
sd.pred <- matrix(NA,14,3)
colnames(sd.pred) <- c("Shift","Lower","Upper")
rownames(sd.pred) <- c("Geographic (CONT)","Geographic (MDN)","Geographic (KNN4)","Geographic (SOI)","Ethnic (MDN)","Ethnic (KNN4)","Ethnic (pSOI)","Trade (MDN)","Trade (KNN4)","Trade (pSOI)","IGO (MDN)","IGO (KNN4)","IGO (pSOI)","Alliance Ties")
for(i in 1:14){
sd.pred[i,]<- sdshift(modz.sd[[i]],sdmean)[2,]
}
sd.pred

## Plot it
pdf("Fig5a.pdf",width=4,height=3)
par(mar=c(2,6,0,0),xpd="n")
plot(x=0,y=0,xlim=c(-.02,.05),ylim=c(0,15),bty="n",axes=F,xlab="",ylab="",col="white")
points(x=sd.pred[,1],y=14:1,pch=3)
#segments(x0=sd.pred[,2],x1=sd.pred[,3],y0=14:1,y1=14:1,col=c(rep("red",4),rep("darkblue",3),rep("darkorange",3),rep("forestgreen",3),"darkgray"),lwd=2) 		# Color
segments(x0=sd.pred[,2],x1=sd.pred[,3],y0=14:1,y1=14:1,col=c(rep("gray30",4),rep("gray65",3),rep("gray30",3),rep("gray65",3),"gray30"),lwd=2) 		# B&W
segments(x0=-.02,x1=.05,y0=0,y1=0)
segments(x0=seq(-.02,.05,by=.005),x1=seq(-.02,.05,by=.005),y0=0,y1=-.5)
segments(x0=0,x1=0,y0=0,y1=15,lty="dotted")
text(x=seq(-.02,.05,by=.01),y=-1,labels=seq(-.02,.05,by=.01),cex=.7)
text(x=.015,y=-2.4,labels="Change in Probability of Transition to Democracy",cex=.7)
text(x=-.05,y=14:1,labels=c("Geographic (CONT)","Geographic (MDN)","Geographic (KNN4)","Geographic (SOI)","Ethnic (MDN)","Ethnic (KNN4)","Ethnic (pSOI)","Trade (MDN)","Trade (KNN4)","Trade (pSOI)","IGO (MDN)","IGO (KNN4)","IGO (pSOI)","Alliance Ties"),pos=4,cex=.7)
dev.off()

###############################
## FIGURE 5b: AIC
###############################

## Extract AIC statistics
aic.mat <- matrix(NA,14,1)
colnames(aic.mat) <- c("AIC")
rownames(aic.mat) <- c("Geographic (CONT)","Geographic (MDN)","Geographic (KNN4)","Geographic (SOI)","Ethnic (MDN)","Ethnic (KNN4)","Ethnic (pSOI)","Trade (MDN)","Trade (KNN4)","Trade (pSOI)","IGO (MDN)","IGO (KNN4)","IGO (pSOI)","Alliance Ties")
for(i in 1:14){
aic.mat[i,]<- modz.sd[[i]]$aic
}
aic.mat

## Plot it
pdf("Fig5b.pdf",width=4,height=3)
par(mar=c(2,.5,0,0),xpd="n")
plot(x=0,y=0,xlim=c(-.02,.02),ylim=c(0,15),bty="n",axes=F,xlab="",ylab="",col="white")
text(x=-.0085,y=-2.4,labels="Akaike Information Criterion (AIC)",cex=.7)
text(x=0.003,y=14:1,labels=round(aic.mat,3),pos=2,cex=.7)
text(x=-.02,y=14:1,labels=c("Geographic (CONT)","Geographic (MDN)","Geographic (KNN4)","Geographic (SOI)","Ethnic (MDN)","Ethnic (KNN4)","Ethnic (pSOI)","Trade (MDN)","Trade (KNN4)","Trade (pSOI)","IGO (MDN)","IGO (KNN4)","IGO (pSOI)","Alliance Ties"),pos=4,cex=.7)
dev.off()

###############################
## FIGURE 5c: Out-of-sample prediction accuracy
###############################

## Create empty matrix to store AUC statistics
aucmat <- matrix(NA,nrow=14,ncol=100)
rownames(aucmat) <- c("Geographic (CONT)","Geographic (MDN)","Geographic (KNN4)","Geographic (SOI)","Ethnic (MDN)","Ethnic (KNN4)","Ethnic (pSOI)","Trade (MDN)","Trade (KNN4)","Trade (pSOI)","IGO (MDN)","IGO (KNN4)","IGO (pSOI)","Alliance Ties")

## AUC LOOP
## * For each of 100 iterations:
## * (a) Randomly select 10% of data for out-of-sample prediction
## * (b) For each of the 14 neighbor definitions:
## *     (i) fit the model on in-sample data 
## *     (ii) calculate the out-of-sample AUC
## *     (iii) save the AUC statistic

for(i in 1:100){
	print(i)
out. <- sample(1:nrow(data.),size=round(.1*nrow(data.)))
out.samp <- data.[out.,]
in.samp <- data.[-out.,]
modz <- models(in.samp)
for(m in 1:length(modz)){
modz.x <- setx(modz[[m]], fn = NULL,data=out.samp)
modz.y <- sim(modz[[m]],modz.x)
modz.fit <- apply(modz.y$qi$ev,2,mean)
aucmat[m,i] <- AUC(out.samp$baut,modz.fit)
}}

## Save output
#save(aucmat,file="AUCMAT.RData")
#load("AUCMAT.RData")

## Calculate mean and 95% CI for 100 iterations of 14 models
auctab <- cbind(apply(aucmat,1,mean),apply(aucmat,1,function(x){quantile(x,.025)}),apply(aucmat,1,function(x){quantile(x,.975)}))

## Plot it
pdf("Fig5c.pdf",width=4,height=3)
par(mar=c(2,6,0,0),xpd="n")
plot(x=0,y=0,xlim=c(.97,1),ylim=c(0,15),bty="n",axes=F,xlab="",ylab="",col="white")
points(x=auctab[,1],y=14:1,pch=3)
# segments(x0=auctab[,2],x1=auctab[,3],y0=14:1,y1=14:1,col=c(rep("red",4),rep("darkblue",3),rep("darkorange",3),rep("forestgreen",3),"darkgray"),lwd=2) # Color
segments(x0=auctab[,2],x1=auctab[,3],y0=14:1,y1=14:1,col=c(rep("gray30",4),rep("gray65",3),rep("gray30",3),rep("gray65",3),"gray30"),lwd=2) # B&W
segments(x0=.97,x1=1,y0=0,y1=0)
segments(x0=seq(.97,1,by=.0025),x1=seq(.97,1,by=.0025),y0=0,y1=-.5)
text(x=seq(.97,1,by=.005),y=-1,labels=seq(.97,1,by=.005),cex=.7)
text(x=.985,y=-2.4,labels="Area Under ROC Curve (AUC)",cex=.7)
text(x=.958,y=14:1,labels=c("Geographic (CONT)","Geographic (MDN)","Geographic (KNN4)","Geographic (SOI)","Ethnic (MDN)","Ethnic (KNN4)","Ethnic (pSOI)","Trade (MDN)","Trade (KNN4)","Trade (pSOI)","IGO (MDN)","IGO (KNN4)","IGO (pSOI)","Alliance Ties"),pos=4,cex=.7)
dev.off()



##########################################################
##########################################################
## Application: The Spread of Democracy                 ##
## Step 3: Simulation of the 2011 Arab Spring           ##
##########################################################
##########################################################



## Run models
modz <- models(data)

## Create list of connectivity matrices
labels <- c("Geography (CONT)","Geography (MDN)","Geography (KNN4)","Geography (SOI)","Ethnic (MDN)","Ethnic (KNN4)","Ethnic (pSOI)","Trade (MDN)","Trade (KNN4)","Trade (pSOI)","IGO (MDN)","IGO (KNN4)","IGO (pSOI)","Alliance Ties")
mats <- list()
mats[[1]] <- mat_cont
mats[[2]] <- mat_dist
mats[[3]] <- mat_knn4
mats[[4]] <- mat_soi
mats[[5]] <- mat_ethmd
mats[[6]] <- mat_ethk4
mats[[7]] <- mat_ethsoi
mats[[8]] <- mat_trmd
mats[[9]] <- mat_trk4
mats[[10]] <- mat_trsoi
mats[[11]] <- mat_igomd
mats[[12]] <- mat_igok4
mats[[13]] <- mat_igosoi
mats[[14]] <- mat_ally

## Run simulations
cnt <- c("Tunisia","Egypt") # Select countries w/ regime change
preds <- regimesim(cnt=cnt,modz=modz,mats=mats,labels=labels)


#############################
## FIGURE 6: Increased risk of regime change (final three)
#############################

cntnum <- which(map$MAP_CNTRY%in%cnt)
MIDs <- map$CNTRY_CODE
ifelse(length(unique(MIDs))==length(MIDs),print("NO DOUBLE COUNTS"),print("DOUBLE COUNTS"))

pdf("Fig6.pdf",height=4.5,width=10) 	# Geographic
layout(matrix(c(1:4,5,5), 3, 2, byrow = TRUE),heights=c(rep(1,2),.25)) 											

## Map loop
nums <- c(1,3,12,14)
for(wnum. in nums){ 							
map$CNTRY_CODE==preds[[wnum.]][,1]
var <- preds[[wnum.]][,"R_AD"]
ci.l <- preds[[wnum.]][,"R_ADl"]
ci.u <- preds[[wnum.]][,"R_ADu"]

## Intervals/colors (color)
#cols <- ifelse(var>=.25,"darkred",
#ifelse(var>=.10 & var<.25,"red3",
#ifelse(var>=.05 & var<.10,"red",
#ifelse(var>=.001 & var<.05,"pink",NA))))
#cols[which(ci.l<=0.001&ci.u>0.001)] <- NA
#cols[map$bautlag == 0] <- NA
#cols2 <- cols
#cols2[-cntnum] <- NA
#cols2[cntnum] <- "black"
#dashes <- 20*(map$baut == 1)

# Intervals/colors (B&W)
cols <- ifelse(var>=.25,"gray25",
ifelse(var>=.10 & var<.25,"gray40",
ifelse(var>=.05 & var<.10,"gray55",
ifelse(var>=.001 & var<.05,"gray75",NA))))
cols[which(ci.l<=0.001&ci.u>0.001)] <- NA
cols[map$bautlag == 0] <- NA
cols2 <- cols
cols2[-cntnum] <- NA
cols2[cntnum] <- "black"
dashes <- 20*(map$baut == 1)

# Plot it
par(mar=rep(.5,4))
plot(map,xlim=c(-7,38),ylim=c(0,45),border="lightblue", col=cols,axes=FALSE)
plot(map, density = dashes, add = TRUE ,col="darkgrey") 
plot(map, border=NA, col=cols2, add = TRUE ) 
legend(x="topleft",adj=c(.2,.4),cex=1.4,col="black",legend=bquote(bold(.(
labels[wnum.]))),box.lty=0)
if(wnum.==1){
	legend(x="topleft",adj=c(.4,.4),cex=1.4,col="black",legend="\n Original",box.lty=0)
	}}

# Draw legend
# fill. <- c("white","darkred","red","red3","darkred")  # Color
fill. <- c("white","gray75","gray55","gray40","gray25") # B&W
labels. <- c("No Change","[.001, .05)","[.05, .10)","[.10, .25)","[.25, 1]")
par(mar=rep(0,4))
plot(0:1,0:1,ylim=c(.8,1),axes=F,col="white")
legend(x="top",cex=1.2,fill=fill.,
bty="n",legend=labels.,
title="Increase in Probability of Democratization",ncol=5)
legend(x="topleft",cex=1.2,fill="black",bty="n",legend="Tunisia and Egypt democratize",title="Counterfactual",xjust=0)
legend(x="topright",cex=1.2,density=c(0,30),
bty="n",legend=c("Democracy","Autocracy"),
title="Initial Regime Type",xjust=.5,ncol=2)
dev.off()


#############################
## TABLE 1: Increased risk of regime change due to Arab Spring
#############################


mat <- preds[[1]][,c("CNTRY_NAME","R_AD")]
mat[preds[[1]][,c("R_ADl")]<0 & preds[[1]][,c("R_ADu")]>0,2] <- 0
for(i in 2:14){
mat[,1+i] <- preds[[i]][,c("R_AD")]
mat[preds[[i]][,c("R_ADl")]<0 & preds[[i]][,c("R_ADu")]>0,1+i] <- 0
}
#mat <- mat[mat$R_AD>0,]
names(mat)[-1] <- labels
mat <- mat[apply(mat[,-1],1,sum)>0.01,]
xtable(mat)
